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The energetics of ionic selectivity in the neuronal sodium channels is studied. A simple model constructed for the 
selectivity filter of the channel is used. The selectivity filter of this channel type contains aspartate (D), glutamate 
(E), lysine (K), and alanine (A) residues (the DEKA locus). We use Grand Canonical Monte Carlo simulations to 
compute equilibrium binding selectivity in the selectivity filter and to obtain various terms of the excess chemical 
potential from a particle insertion procedure based on Widom's method. We show that ions in competition 
with Na”^ are efficiently excluded from the selectivity filter due to entropic hard sphere exclusion. The dielectric 
constant of protein has no effect on this selectivity. Ca^”^ ions, on the other hand, are excluded from the filter 
due to a free energetic penalty which is enhanced by the low dielectric constant of protein. 
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l. Introduction 

Sodium (Na) channels can be categorized on the basis of their function, the cell in which they are 
found, structure of the protein (both secondary and tertiary), and the structure of the selectivity filter 
(SF). The SF is a narrow region of the permeation pathway, where the channel discriminates between 
different ions. The selectivity properties of different channels primarily depend on the type of amino 
acid motifs present in their SF. 

The two most widely studied classes of Na channels are the neuronal (this is the one studied here) and 
bacterial Na channels. The SF of neuronal Na channels has a DEKA locus made of aspartate (D), glutamate 
(E), lysine (K), and alanine (A) residues. On the basis of their homology with L-type calcium (Ca) channels 

m, these amino acids seem to face the permeation pathway. The accurate structure of the DEKA Na 
channels is still unknown, so theoretical studies are restricted to using models based on homologies of 
the known structures or on reduced models based on minimal structural information available. This is 
the approach used in this work, while the minimal structural information is that the SF has the DEKA 
locus. 

The bacterial Na channels, on the other hand, have X-ray structures measured recently (2}{5l. These 
channels include NavMs (2] El, NavAb (3]|4l, and NaChBac m The structure of a Ca^”^-selective mutant 
of NavAb is also available m These channels have a lot of aspartates and glutamates in their SF. There¬ 
fore, at a first glance, they look like a Ca channel. Hydration plays an important role in the selectivity 
mechanisms of these channels (81, but this is not the subject of the present study. 
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Simulation studies for Na channels have been based on models of different resolutions. All-atom, 
explicit-water models are usually used when X-ray structures are available. They are generally studied 
using molecular dynamics (MD) simulations |1|[5]|9). In the case of the DEKA channel, Lipkind and Foz- 
zard ITOl performed MD simulations to explore the Na”^ vs. selectivity for various mutants of DEKA 
based on extreme homology modelling. 

Boda et al. 11111131 and Vora et al. QH used reduced models of Na channels in the implicit solvent 
framework. In these models, only the SF amino acids were represented in an explicit way, while other 
parts of the channel protein were reduced into a dielectric body. This is the modelling level that we 
use in this work. An intermediate approach is that of Finnerty et al. (15), who proposed a localization 
method, where SF aminoacid terminal groups are localized into certain positions inspired by structural 
information. 

The advantage of reduced models is that they permit the design of simulation setups in time and 
length scales that mimic experimental setups and are capable of studying a wide range of concentrations 
and voltage. Also, they make it possible to focus on the essential features of the system (SF structure, pore 
geometry, bath concentrations, voltage, etc.) and to take the effect of the remaining degrees of freedom 
into account in an averaged, but physically well-based manner (dielectric response as well as external 
constraints such as the walls of the channel and the membrane). 

Simulations can also be distinguished on the basis of the fact whether they were performed in or out 
of equilibrium. Equilibrium simulations can study the selective binding of various ions in the SF. Monte 
Carlo (MC) simulations, especially in the grand canonical (GC) ensemble (Grand Canonical Monte Carlo, 
GCMC) are ideal tools for this purpose 11111131 [T^ . MD simulations can also be used to study selective 
binding (To][T6][T3. Some properties of transport, however, can be extrapolated even from equilibrium 
simulations on the basis of the integrated Nernst-Planck equation as suggested by Gillespie et al. (18] 
[19). Simulating transport requires a dynamical simulation method. These can be MD simulations (1][5] 
|9), Brownian Dynamics simulations (T4]|20H22), and Dynamical Monte Carlo (DMC) simulations (T^l^ 
El. Transport can also be studied with theoretical methods such as the Energy Variational approach of 
Eisenberg et al. 12511281 . 

To extend equilibrium binding-selectivity simulations to non-equilibrium situations of steady-state 
ionic transport is of crucial importance because experimental data are available for ionic currents from 
electrophysiological measurements (291144) . The relation of the fluxes carried by the competing ions (flux 
ratio) deflnes dynamical selectivity. How binding selectivity is related to dynamical selectivity is, how¬ 
ever, a non-trivial issue as shown by Rutkai et al. (H. In particular, the flux is determined not only by 
the occupancy of a given ionic species in the channel, but also by its mobility. 

Measurements show permeability ratios 0.06 and 0.13 for K'^/Na'^ and Ca^'^/Na'^, respectively f29\ - 
(38), while < 0.01 flux ratio for (39ll44l . To a first approximation, we can assume that binding 

selectivity agrees well with the above selectivity values measured in terms of flux. To what degree this 
assumption is valid can be studied by dynamical simulation methods. Our first attempt in this direction 
is the DMC study of Csanyi et al. (13). 

In this paper, we focus on equilibrium binding, so it is a direct continuation of our previous papers 
(TTI (12), where the binding selectivity of the DEKA locus was studied with GCMC simulations using a 
reduced model of the SF. These studies used the charge-space competition (CSC) mechanism of Nonner 
and Eisenberg 14611521 extended later to inhomogeneous models of the channels studied by GCMC simu¬ 
lations (TT1IT21IT8IIT9II53H58). 

The main conclusions of those papers (TT][T2) were that K”^ ions are excluded from the SF by steric 
repulsion, while Ca^”^ ions are excluded by an electrostatic penalty. The new aspect of this study is that 
we provide an energetic analysis for the phenomena described in our 2007 paper (12). The energetic 
analysis is performed by separating the free energy (more exactly, the chemical potential) into vari¬ 
ous terms corresponding to various interactions such as volume exclusion, ion-ion, ion-dielectrics, self¬ 
energy, etc., interactions. This approach was introduced by Gillespie (59) (60) in his density functional 
studies for the Ryanodine Receptor Ca channel and was extended to three-dimensional models including 
inhomogeneous dielectrics using a GCMC methodology (61) on the basis of Widom’s particle insertion 
method (62]|63. 

In our previous work, we analyzed the energetics of the selectivity of the F-type Ca channel (61). In 
that paper, the dielectric constant of the protein (Cpr) was allowed to be different from that of the baths 
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(Cw). It was shown that the low dielectric protein surrounding the pore focusing the electric field, and 
thus enhancing the electrostatics, is necessary to reproduce the micromolar selectivity observed for the 
L-type Ca channels (64l|65). We also extended that work for the case of a dielectric constant different 
inside the channel (Cch) from that of the bath This model is a simple representation of solvation. Our 
results showed that solvation plays a minor role in the selectivity mechanism of the L-type Ca channel. 
The explanation is that the solvation penalty for Ca^”^ is balanced by stronger interactions of Ca^”^ with the 
SF charges. Our simulations extending this work to the DEKA locus are in progress and will be published 
in a subsequent paper. 

In this work, however, we restrict ourselves to the case, where the dielectric constants inside and out¬ 
side the channel are the same (Cch - This is the model that was considered in our work (12) in 2007. 
The SF of the Ca channel is highly charged (EEEE locus, four glutamate residues providing -Ae charge). 
The DEKA filter, on the other hand, is weakly charged (-le altogether). Therefore, it does not favor diva¬ 
lent ions (Ca^”^). Additionally, the bulky terminal group of the lysine is present, which, according to our 
hypothesis, is there to exclude large ions such as This paper examines how these mechanisms work 
and their energetic basis. 


2. Model 

In our model, most of the atomic structure of the Na channel is reduced to a coarse-grained geometry 
(figure]^. The channel protein is represented as a continuum solid with dielectric coefficient Cpr. The 
three dimensional body of the protein is obtained by rotating the thick line in figure about the r = 0 
axis. The protein thus forms an aqueous pore that connects the two baths. Water in the baths and pore 
is described as an implicit solvent that is a continuum dielectric having uniform dielectric coefficient 
Cw = 80. The central cylindrical part of the pore (with radius i? = 3.5 - 4.5 A and length 10 A) forms the 
selectivity filter that includes the only atoms of the protein that are treated explicitly. These atoms are 
four half-charged 'oxygen ions’ [figure[^(B), red spheres] representing the charged terminal groups 
of the D and E residues, while a positive 'ammonium ion’ NH^ [figure [^(B), blue sphere] represents the 
terminal group of the K residue. The alanine is ignored. The structural oxygen ions are confined to the 
selectivity filter (their centers are in the region r ^R-Rf, |z| ^ SA-Rf, where Rf is the ionic radius), but 
they can move freely inside the filter. 

The ions are modelled as charged hard spheres with crystal radii (see caption of figure [^. The com¬ 
putation of the intermolecular energy terms due to screened Coulomb potentials and interactions with 
polarization charges induced on the dielectric boundaries [the boundary of the protein and the elec¬ 
trolyte; thick line in figure (A)] are described in our previous works (55] [611 [63. Ions are restricted to 
the aqueous space of the model and cannot overlap with hard walls in the system. Figure (A) shows 
only the small central region of the simulation cell. The entire simulation cell is a cylinder having typ¬ 
ical dimensions of radius 40 A and length 180 A. The channel is embedded in a membrane region that 
excludes ions by hard walls as described before (55). 


3. Method of energetic analysis 

In an equilibrium GCMC simulation, the acceptance of ion insertion/deletions of ions is governed by 
the configurational chemical potential of the respective ionic species i defined as 

Pi = /crinQ(r) + juf^(r) = /cTlncKB) -E juf^(B), (3.1) 

where k is Boltzmann’s constant, T is the temperature, Ci(r) is the concentration profile, jU?^(r) is the 
excess chemical potential profile, Ci(B) is the bulk concentration, and jU?^(B) is the bulk excess chemical 
potential. Although kTci (r) and jU?^(r) can be different in different regions (they are position dependent), 
their sum is constant due to equilibrium. The bulk excess chemical potentials corresponding to the 

prescribed bulk concentrations Ci(B) are calculated using the Adaptive GCMC method (^. By rewriting 


13601-3 




D. Boda etal. 


15 


10 


5 - 


-5 - 


-10 


-15 




-15 -10 


n ^ r 


protein pr 


protein 


J_^_L 


-5 0 

z / A 




10 15 



Figure 1. (Color online) Model of ion channel, membrane, and electrolyte. The three-dimensional geome¬ 
try (B) is obtained by rotating the two-dimensional shape shown in panel A around the z-axis. The simula¬ 
tion cell is much larger than shown in the figure. The blue lines represent the grid over which the excess 
chemical potential profiles are computed. The grid is finer inside the channel (width 0.5 A), while it is 
coarser outside the channel (width 2 A). The selectivity filter (|z| < 5 A) contains 4 half charged oxygen 
ions (red spheres in panel B) and an ammonium ion NH^ (blue sphere in panel B). For the radii of 
the ions, the Pauling radii are used: 0.6, 0.95, 1.33, 1.52, 1.7, 0.99, 1.81, 1.4, and 1.5 A for Li”^, Na”^, 
Rb”^, Cs”^, Ca^”^, CP, and NH^ respectively. 


equation l |3.1) , the excess chemical potential difference is defined as 


(r) = W - (B) = -fcrin 


Ci(r] \ 

Ci (B)j' 


(3.2) 


It can be identified with the binding free energy of an ion moved from a bath (B) to position r of the 
channel (61). If we write up equation l |3.2| for Na”^ and and take the difference, we can derive that 


( W ] _ f CNa^ (B) ] ^ 
icK+(r)j '^Uk+CB)] kT 


(3.3) 


where 

AA;U^^(r) = AMNa+ W - A/u^ (r). (3.4) 

Similar equations can be given for other pairs of ions. 

In equation (3.3) , the left-hand side is called 'binding selectivity’ because it expresses the degree to 
which Na”^ is favored over at location r (binding selectivity is positive if location r is selective for Na”^ 
over K”^). The corresponding term on the right-hand side containing the bulk concentrations is called 
'number advantage’ (59) because it expresses the advantage that an ionic species gets from outnumbering 
the other species in the bulk. The channel can become selective for a given ionic species for two reasons: 
either from the number advantage or the energetic advantage expressed by AAju^^(r). 

The energetic advantage, however, contains terms due to different interactions present in the system 
as described in appendix]^ The EX term can be divided in various ways. Here, we use the division used 
in our latest work (66): 

AjufX(r) = Ajuf^ (r) -f Ajuf (r) -f Aju]^ (r) -f (r) (3.5) 
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or briefly EX = HS + II + ID + SELF, where HS means hard sphere exclusion, II means interaction with the 
ions, ID means interactions with the dielectrics (polarization charges induced by other ions), and SELF 
means interactions with the polarization charges induced by the ion itself. (In the division of our earlier 
work 1611. we used the DIEL term that contained the SELF term, namely, DIEL = ID + SELF.) We can also 
use the division EX = HS+MF+SC+SELF, where MF means the interaction with the mean (average) elec¬ 
tric field of all the existing charges in the system (ionic and induced). SC expresses correlations beyond 
the mean field level (SC refers to 'screening’) (^. The SELF term is a one-particle term (mean-field in 
nature) and corresponds to the average electrostatic interaction energy of the inserted ion with its self- 
induced charge. It is not included in the ID or the MF term. The SELF term corresponds to the dielectric 
boundary force or energy of reference (6^ . 

The computation of all these terms can be found in our original paper |6T) and in appendix]^ Briefly, 
the total EX chemical potential can unambiguously be obtained by inserting charged hard spheres (repre¬ 
senting the ions) in the Widom particle insertion method. Different terms of EX are computed by inserting 
particles interacting only through short-ranged (HS) or more direct (II) interactions and obtaining the rest 
as residuals. For example, it is reasonable to compute the HS term by inserting uncharged hard spheres 
with the same radius as the respective ion in the Widom procedure. All the remaining terms (II, ID, SELF) 
are electrostatic in nature and obtained by deducting the HS term from the EX term. (The separation of 
HS and electrostatic terms and their effect on selectivity can already be found in the work of Nonner et 
al. in the context of the mean spherical approximation.) Similar procedures are used to separate the 
II and ID, as well as the MF and SC terms, as described in appendix]^ 

The r-dependence of various terms is computed by ion insertions into grid cells shown in figure 
Note that the concentration profile can be computed in two different ways. First, sampling the number 
of ions in a volume element, computing the average ion number and dividing by the volume of the ele¬ 
ment. This is advantageous when the concentration and/or the volume element is large so there is a large 
enough sample of ions. The concentration, on the other hand, can also be computed from equation j3.1| l 
by computing the EX term from the Widom method and deducting it from the chemical potential. This 
approach is useful where the concentration is low. This method was used in our simulations for the DEKA 
channel. 

Our grid is two-dimensional because we have rotational symmetry. Our profiles, therefore, are ex¬ 
pressed in terms of the (z, r) cylindrical coordinates. In this work, however, we show the results that are 
averaged over the r-coordinate 



(3.6) 


where Rminiz) - Riz)- i?iargerion(-^) is the cross-section that is accessible to the center of the larger of the 
competing ions, [Riz) denotes the radius of the simulation domain at z]. 

4. Results and discussion 

We start our discussion with competition of ions of the same charge. Specifically, we study the selec¬ 
tivity of Na”^ over various monovalent ions. In the classical mole fraction experiment, the mole fraction 
of one ion (Na”^, for example) is changed while keeping the total cation concentration constant (when 
divalent ion is present, the total ionic strength is kept constant in some studies). These results are seen in 
figure 5 of reference flTi . In this work, the concentration of the two competing cations in the baths is the 
same (50 mM), so the number advantage is zero. 

Figure shows the various terms of the A/i?^(z)-profiles for Na”^ and for protein dielectric con¬ 
stant Cpr = 10 and filter radius R = 3.5 A. The value Cpr = 10 is the value fixed in our studies for the L-type 
Ca channel (T8][T9][55l[56]|58][6Tl[70l|7^ The value i? = 3.5 A was used in our DMC study for the DEKA Na 
channel to reproduce experimental data 1131 . 

The EX terms are related to the concentration ratios through -ln[c/(r)/c/(B)] [see equation 13.21 1. 
Therefore, where the EX term (or any component) is negative, it energetically favors the ionic species, so 
it increases the concentration of that ionic species. As also seen in figure 6 of our previous paper |T2l, 
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50 mM Na"" vs. 50 mM K"" 



Figure 2. (Color online) The Aju^-(z)-profiles for Na”^ and for the case when the bath concentration is 
the same for the two competing monovalent cations (50 mM) for Cpr = 10 and R = 3.5 A. 


there are peaks at the entrances of the SF and the vestibules (,\z\ ~ 5 A). In the center of the SF, on the 
other hand, the concentrations are low. This region forms a depletion zone for both ions, where ions have 
difficulty to enter. The question, therefore, is which ion is excluded less from this region. The answer is 
that there are more Na”^ than in the SF (the EX term is lower for Na”^), so the SF is Na'^-selective. 

All the electrostatic terms (II, ID, MF, SC) are negative except the SELF term. The SELF term is repulsive 
because the ions are in the = 80 region, so the sign of the induced charge on the CprlCw boundary is 
the same as the sign of the inserted ion itself. This practically corresponds to the dielectric penalty an ion 
must pay when it passes the low dielectric membrane region as described in classical works 17211751 . The 
SELF term is slightly larger for Na”^ because the smaller Na”^ can get closer to the channel wall and can 
induce a larger polarization charge. 

The other term that is positive is the HS term describing the volume exclusion. This is the term that 
is very different in the case of Na”^ and it is larger in the case of Since the size of ions (we 
talk about the dehydrated (Pauling) radius) is larger, it is more difficult to insert such an ion in the SF. 
Therefore, has a larger entropic penalty than Na”^ does. This difference is especially apparent in the 
center of the SF, where the NH^ (the structural ion representing the large terminal group of the lysine) 
profile has a peak (see figure 6 of reference (l2l). Without the HS term (ions of finite size) we could not 
get a Na”^-selective filter (against K”^) in this model. 

The MF term is negative, because the SF is negatively charged. There is no space for the cations to 
fully neutralize the SF charge. The SC term is similar to the MF term in order of magnitude indicating 
that mean field theories are not sufficient to study ionic systems in crowded confined spaces such as the 
SF of ion channels. 

The dominant term that drives Na”^ vs. selectivity is the HS term. In figures [s] and therefore, 
only the differences of the EX and HS terms are shown for various cases. In this special case, where the 
number advantage is zero, the EX difference is equal to the binding affinity [see equation l |3.3) ], while the 
HS term is the dominant term of AAjU^^(r). Since the differences are obtained by deducting the terms 
from the Na”^ terms, positive values favor Na”^. 

Figurej^shows the profiles for various pore radii for Cpr = 10 (top panels) and Cpr = 80 (bottom panels). 
Narrower channels favor Na”^ even more, as expected, because it is even more difficult to find space for 
the large ions in the small SF compared to Na”^. Putting it in another way, Na”^ vs. selectivity is 
better for narrow channels, where stronger competition is forced by the confinement and lack of space. 
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Binding affinity HS advantage 



Figure 3. (Color online) The binding affinity and HS advantage curves for Na”^ vs. competition for 
epr = 10 (top panels) and 80 (bottom panels) for filter radii R = 3.5, 4, and 4.5 A (50-50 mM bath concen¬ 
trations). 


Binding affinity HS advantage 



z / A z / A 


Figure 4. (Color online) The binding affinity and HS advantage curves for the competition of Na”^ against 
various monovalent ions (Li”^, Rb”^, Cs”^) for Cpr = 10 (top panels) and 80 (bottom panels) (50-50 mM 

bath concentrations, R = 3.5 A). 
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so the smaller size of Na”^ has the advantage. The binding affinity curves (left-hand panels) and the HS 
advantages (right-hand panels) behave similarly with small differences due to other energetic terms (see 
figure [^. 

Another conclusion of the figure is that Na”^ vs. selectivity does not depend on the dielectric con¬ 
stant of the protein; the curves for Cpr = 10 (top panels) and Cpr = 80 (bottom panels) behave practically 
the same. 

Figure|^shows the same curves but now for a fixed pore radius (R = 3.5 A) and different monovalent 
cations (Li”^, Rb”^, Cs”^) competing with Na”^. The main conclusion is similar to those drawn in figurej^ 

the crowded SF favors the smaller ion. The pore is selective for Li”^ against Na”^, while it is selective for 
Na”^ against the larger ions. 

The protein dielectric constant does not have an effect on these profiles. Of course, the value of Cpr 
has a large effect on the individual ionic profiles and the occupancies (see figure 8 of reference (12)), but 
not on the relative ones that we study here. 

In the second half of this section, we analyze the competition of Na”^ against Ca^”^. The other usual way 
to study the behavior of the channel having varying electrolyte composition is to keep the concentration 
of one species fixed (Na”^, for example) and to add another species (Ca^”^, for example) gradually. This 
added salt experiment was done by Aimers and McCleskey in their experiment for the L-type Ca channel 
(64l|65). We performed this kind of experiment in our previous simulations for the DEKA locus and its 
DEEA mutant, see figure 2 of reference IT^ . 

Those simulations qualitatively reproduced the experiment of Heinemann et al. (38). Heinemann 
et al. found that mutating the DEKA locus into a DEEA locus, the selectivity behavior of the channel 
is reminiscent to Ca channels rather than Na channels. In experiment, the current drops to half (IC50) 
at Ca^”^ concentration 10“^ M, while in our simulations, the number of Na”^ ions drops to half at the 
same concentration. The explanation is that the DEEA mutation has -3e charge producing a Ca channel, 
but with weaker selectivity than in the case of the -Ae charge (EEEE locus). The DEKA locus, on the 



Figure 5. (Color online) Ca^”^ and Na”^ concentration profiles for two different Ca^”^ concentrations (10 
and 40 mM in top and bottom panels, respectively) with a 50 mM Na”^ background (R = 3.5 A). The profiles 
are shown for protein dielectric constants Cpr = 80 (left-hand panels) and 10 (right-hand panels). 
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other hand, shows over Ca^”^ selectivity. This selectivity is stronger for smaller Cpr (see figure 10 (A) 
of reference El). The dielectric constant of the protein, therefore, has a strong effect in the case of 
monovalent vs. divalent competition. 

In figure we show the results only for two chosen Ca^”^ concentrations, 10 mM (top panels) and 
40 mM (bottom panels) — both are well above the physiological values (~l-2 mM). 

The background Na”^ concentration is 50 mM. The Na”^ and Ca^”^ concentration profiles are shown for 
Cpr = 80 (left-hand panels) and 10 (right-hand panels). 

There are more Na”^ than ions in the filter in the case of Cpr = 10 for both concentrations. A single 
Na”^ ion efficiently counterbalances the filter charge. Ca^”^ ions, on the other hand, overcharge the filter, 
which is electrostatically unfavorable. To counterbalance this overcharge, a cr would be needed, but 
there is no space left for it in the filter. 

In the case Cpr = 80, on the other hand, there are more Ca^”^ ions at [Ca^”^] = 40 mM. The explanation 
is that Ca^”^ is still double charged so the SF attracts it more strongly. The overcharged filter is balanced 
by cr ions from outside the filter. In this case, this is possible because the Coulomb forces are more 
long-ranged and more screened than in the case of Cpr = 10, where the low-dielectric protein focuses the 
electric field. This means that the low dielectric protein is needed to exclude Ca^”^. 

The energetics of this phenomenon is analyzed in figures and The difference in Na”^ vs. Ca^”^ 
selectivity is more clearly seen by plotting the binding selectivity curves. When this is positive, the pore 
is Na”^-selective, while it is Ca^”^-selective in the opposite case. The number advantages are also indicated 
with dashed horizontal lines. As bath Ca^”^ concentration is increased, this line and the binding selectivity 
curve with it are shifted downwards. The shape of the binding selectivity curves does not change much 
with the bath Ca^”^ concentration. We can conclude, therefore, that Na”^ vs. Ca^”^ selectivity does not 
depend on the bath Ca^”^ concentration. This is because the DEKA locus is a singly occupied SF; only one 
cation occupies the SF at one time (or none). 

This was not true for the L-type Ca channel. That channel could be multiply occupied, so selectivity 
behavior was a function of Ca^”^ concentration due to correlations of cations in the filter. Furthermore, 
the SF of the EEEE locus became more charge neutral as Ca^”^ concentration was increased. Because of 
that, the MF terms decreased (see figure 7 of Boda et al. (gT)). That effect is absent here; the probability 
that a channel becomes charge neutral does not depend on ionic concentrations, but it rather depends on 
entropic effects (available space in the channel given by filter radius and ion sizes). 


Binding selectivity, added to 50 mM Na^ 



Figure 6. (Color online) Binding selectivity (solid lines) and number advantage (dashed lines) curves for 
the cases considered in figure[^for three different concentrations. 
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Figure 7. (Color online) The various terms of the free energy advantage for [Ca^'^]=10 mM and 
[Na'^]=50 mM. Top and bottom panels show the results for epr = 10 and 80, respectively. Left panels 
show the EX curve in the EX = HS + II -h ID -h SELF division, while the right panels show the EX curve in 
the EX = HS + ME -h SC + SELF division. 


The difference of bi ndin g selectivity and number advantage defines the free energy advantage, 
AAjLi^^(r), [see equation The terms of that advantage are analyzed for Cpr = 10 and 80 for a given 
Ca^”^ concentration (10 mM) in figurej^ 

The top panels show the Cpr = 10 data. The left-hand panel shows the II and ID terms (EX = HS + II -i- 
ID -f SELF), while the right-hand panel shows the MF and SC terms (EX = HS + ME -f SC + SELF). The EX 
and SELF terms are shown both in the left-hand and right-hand sides. The HS term (not shown) is close 
to zero because the ions are of similar sizes. The EX term is also close to zero in this case, but this is the 
effect of the balance of the different free energy advantage terms. The SELF term is very positive, so it 
favors Na”^. This term is about four times larger for Ca^”^ than for Na”^ so it plays the role of solvation 
penalty in this model. Without the SELF term, we could not get a Na”^ selective filter (against Ca^”^) in this 
model. Both the II and ID terms (as well as the MF and SC terms, see right-hand panel) favor Ca^”^ because 
Ca^”^ is attracted twice as strongly by the SF charges (ionic and induced) as Na”^. 

The bottom panels show the Cpr = 80 data. Here, the ID and SELF terms are absent, because there 
is no dielectric boundary present. The ID term favors Ca^”^, while the SELF term favors Na”^. Since the 
SELF term is larger in absolute value, these two terms together (ID -i- SELF) still favor Na”^, so the channel 
becomes less Na”^ selective in their absence. 

The SC term is small for Cpr = 80, which means that Na”^ vs. Ca^”^ selectivity is chiefly a mean-field 
effect in this case; the 0^^^“ ions attract Ca^”^ twice as strongly as they attract Na”^. In the case of Cpr = 10, 
on the other hand, SC is quite large indicating a SF of higher density and correlations beyond the mean- 
field level (mainly, with induced charges). 

Summarizing, the EX term is negative for Cpr = 80, so it is rather a Ca channel. The EX term is close 
to zero for Cpr = 10, which means that neither ions are favored energetically. Binding selectivity is driven 
by the number advantage, which results in a Na”^ selective channel (against Ca^”^) at physiological Ca^”^ 
concentrations (1-2 mM). 
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5. Conclusions 

We analyzed the energetics of ion selectivity in the SF of the DEKA Na channels. The reduced model 
studied before (13 was capable of reproducing the basic characteristics of this channel. We showed that 
ions are excluded from the SF due to entropic hard sphere exclusion. The dielectric constant of the 
protein has no effect on this selectivity. In general, this filter favors smaller ions over larger ones. 

Ca^”^ ions, on the other hand, are excluded from the filter due to a free-energetic penalty which is 
enhanced by the low dielectric constant of the protein. The DEKA locus works as a Na channel in the 
Na”^ vs. Ca^”^ competition by not favoring Ca^”^. The dominant term is the number advantage in the bulk 
solutions. In physiological situations this mechanism suffices. 

We showed that the dominant term of the energetic penalty is the SELF term, which is a dielectric 
penalty — the interaction of the ion with the polarization charges induced by itself. This dielectric penalty 
is a simple, implicit representation of solvation penalty in the framework of this model, where Cch = ^w- 
Simulations, where a different dielectric constant inside the channel is used, take solvation into account 
explicitly. 
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A. Widom particle Insertion method to compute the components of the 
excess chemical potential 


The excess chemical potential profile can be computed using Widom’s particle insertion method 
(^[ 63 . We divide the simulation cell into small volume elements as described in reference (61) and 
insert “ghost” particles into uniformly generated positions in these volume elements. We compute the 
interaction energy U (r) of the “ghost” ion inserted at position r with the whole system and use it in the 
operation 

r [[/(r)] (A.l) 

where the brackets denote GC ensemble average. If the interaction energy U (r) contains all the terms (no 
matter how it is divided into terms), operator W provides the full excess chemical potential 

^^(r) = W [J7P(r) + J7!’(r) + uf{r) + . (A.2) 

A diverging term corresponding to overlap with protein and membrane walls is omitted in this 

equation, because we evaluate the excess chemical potential only at the allowed positions. 

The II term of the energy is obtained as Ufir) = ZiZje^y/^^jir^rj), where 


Vliri.rj] ■■ 


BTrcoCwlr^-r/ 


(A.3) 


describes the Coulomb interaction between two unit charges at positions r/ and rj. The ID term is ob¬ 


tained as uP{r) = ZiZjey/^^ir,rj), where 


IDr 




Bttco 


J |r/-s| J |r;-s| 


(A.4) 
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describes the interaction of a unit charge at r/ with the polarization charge, hjirjyS), induced by another 
unit charge at rj (or vice versa). Vector s runs over the dielectric boundary The polarization charge is 
determined using our Induced Charge Computation method I55II671 . 

We define the terms in the excess chemical potential that correspond to different interactions as sug¬ 
gested by Gillespie (^. The definition of these terms is not unique. In our previous work (61), we sug¬ 
gested a possible and physically well-based procedure. The HS term in the excess chemical potential is 
computed by inserting uncharged hard spheres into the system with the same size as the corresponding 
ion, but without the charge: 

nf^{r)^W[uf^{r)]. (A.5) 

The II -I- ID -I- SELF part is the difference EX - HS. If we insert charged hard spheres into the system, but 
ignore their interactions with the polarization charges, we can compute an excess chemical potential 
term describing the ion-ion interactions including the HS interactions: W [[/P^(r) -F t/!^(r)]. The II term 
(that corresponds solely to the interaction with the ionic charges) is obtained by subtracting the HS term: 

l/l (r) = r [ (r) + (r) ] - r [ (r) ]. (A.6) 


The ID term (that corresponds to the interactions with polarization charges induced by other ions) is what 
remains: 


M!°(r) = A<r(r)- 




(A.7) 


The SELF term is a one-particle term that corresponds to the i = j term of the ID energy in equation l |A.4| . 
The MF term is simply the interaction with the mean electric field computed by sampling with a unit point 
charge as described in reference (61). The SC term, again, is what remains: SC = EX - HS - MF - SELF. 
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EHepreTi/iKd iohhoY KOHKypeHL^iY y ce/ieKTi/iBHOMy (|)mbTpi DEKA 
HeMpOHHI/IX HdTpiCBI/IX KdHd/IIB 

BoflaR r. J]\(SP, Ax. OoHceKcR B. AM3eH6epi® 

^ OaKy/ibier (t)i3i/iMHOi xiMii, yHiBepci/iTer riaHHOHii, BecnpeM, H-8201, yropu^i/iHa 
^ BiA^meHHfi MaTeMaiMKi/i i KOMn'HDiepHi/ix Hayx, AproHHCbKa Hai4ioHa;ibHa yiaOoparopm, 

M. AproHH, IL 60439, CLIJA 

^ Mepexa o6Mi/ic/iK)Ba/ibHoT HaHOTexHO/iorii, yHiBepci/iTer flepAbK), m. BecT-/la(t)aMeTT, mraT iH^iaHa, CLUA 
^ OaKy/ibier MO/iexy/i^ipHoT 6io(t)i3i/iKi/i i (t)i3io/ioriT, MeAi/iHHi/iM LteHip yHiBepci/iTery Pama, 

By/i. Beer fappicoH, 1750, Mixaro, CLUA 


ripoBeAeHo Aoc/iiAxeHHJi eHepren/iKi/i ioHHoi ce/ieKTUBHocri b HewpoHHi/ix xaHa/iax HaipiK). Bi/iKopi/iciaHo npo- 
ciy Mo^e/ib, cKOHcrpywoBaHy cnepiayibHo ce/ieKTUBHoro (t)i/ibTpy. Ce/ieKii/iBHi/iM (|)mbTp xaHa/ibHoro n/iny 
MicTi/iTb 3a/ii/iLiJKi/i acnapary (D), r/iKuaMary (E), /liai/iHy K) ra a/iaHiny (A) (o6/iacTb DEKA). Bi/iKopi/icraHO Mo^e- 
/iK)BaHH?i MeroAOM MoHie Kapyio y Be;ii/iKOMy xaHOHiMHOMy aHcaM6;ii o6Hi/ic/ieHHJi ceyieKTUBHocii piBHo- 
BaxHoro 3B'?i3yBaHH?i y ce/ieKiuBHOMy (|)i/ibTpi i oipuMaHH^i pi3Hi/ix H/ieHiB HaA/imuKOBoro xiMiMHoro no- 
TeHi4ia/iy b peay/ibiari npoLteAypi/i BCTaB/ifiHH?i MacTHHOK Ha ochobI Mero^y Bi^oivia. rioKaaaHo, uxo ioHi/i y 
cynepHHLtTBi 3 Na”^ e(t)eKTi/iBHo Bi/i/iyMaK)Tbc?i 3 ce/iexTi/iBHoro (t)i/ibTpa 3a paxyHOK eHiponiMHoro Bi/iK/iK)HeHHJi 
TBepAHX c(t)ep. /l,ie/ieKTpi/iMHa npoHi/iKHicib npoTemy He Mae xoAHoro BnyiHBy na Aany ceyieKTHBHicib. 3 inujoro 
6oKy, ioHi/i Ca^”^ Bi/i/iyMaK)TbCfi 3 (t)i/ibTpa 3a paxynoK Bi/ibHoro enepreTi/iMHoro "nena/ibii", uxo niACi/i/iweTbca 
HH3bKOK) Aie/ieKipHHHOK) npoHHKHicTK) npoieiHy. 

KniOHOBi c/ioea: MoHre Kap/io, npi/iMiri/iBHa MOfl,enh eneKiponiriB, ioHHviiA Kanan, ce/ieKTi/iBHicTb 
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